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Preface 


'•  T  The  purpose  of  this  study  was  to  investigate  the  effect 

a  trailing  vortex  wake  has  on  an  airfoil  undergoing  a  con- 

AiS  y  pV\  y  5 

>  in  two-dimen¬ 
sional,  incompressible,  irrotational  flow.  Potential  flow 
theory,  conformal  mapping  by  the  Joukowski  transformation, 
and  numerical  integration  and  differentiation  techniques  were 
used  to  develop  a  computer  algorithm  to  model  the  problem. 
Once  the  program  was  formulated,  it  was  used  to  solve  the 
impulsive-start  problem  of  airfoil  motion.  The  results  were 
found  to  be  in  excellent  agreement  with  the  results  obtained 
by  others.  When  applied  to  the  constant  rate-of -change  of 
angle-of -attack  problem,  the  results  showed  that  a  trailing 
vortex  wake  has  a  measurable  and  predictable  effect  on  the 
production  of  lift  on  an  airfoil  undergoing  a  constant  £\>J  ■ 

iiu,. 

^  /Che  results  of  this  work,  taken  alone,  are  helpful  in 
understanding  the  phenomena  known  as  dynamic  stall,  tout^  ’  • 

coupled  with  existing  boundary-layer  studies  the  results 
may  lead  to  additional  understanding  of  the  phenomena.  More 
specifically,  the  computer  program  developed  here  could  be 
used  to 'more  realistically  predict^the  inviscid  flow  about 
a  pitching  airfoil  as  it  approaches  the  dynamic-stall  con¬ 
ditions  . 

This  study  could  never  have  been  completed  without  the 
help  of  others.  I  owe  a  great  deal  of  thanks  to  Major  Eric 
Jumper,  who  not  only  posed  this  problem  to  me,  but  who  also 


provided  invaluable  assistance  throughout  the  investigation. 
I  also  wish  to  thank  Lt .  Colonel  Michael  Smith  for  his  help 
with  potential  flow  theory  and  his  expert  advise  regarding 
the  writing  of  this  report.  Finally,  I  wish  to  thank  my 
wife  Anneliese  for  her  translations  of  the  German  references 
and,  most  importantly,  for  her  help  in  keeping  the  work  in¬ 
volved  in  this  thesis  in  proper  perspective. 


Kenneth  W.  Tupper 
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Abstract 

This  study  explored  the  effect  of  a  trailing  vortex 
wake  on  the  production  of  lift  on  an  airfoil  undergoing  a 
constant  rate  of  change  of  angle  of  attack,  a.  The  study 

showed  that  when  an  airfoil  encounters  a  constant-a  flow, 

/ 

the  trailing  vortex  wake  acts  to  suppress  the  slope  of  the 
airfoil's  vs.  a  curve.  The  change  in  magnitude  of  this 
effect  as  a  function  of  airfoil  thickness  and  camber  was 
also  investigated. 

Potential  flow  theory  was  used  to  model  the  flow  about 
a  two-dimensional  circular  cylinder,  and  that  flow  was  trans¬ 
formed  to  flow  about  an  airfoil  by  the  Joukowski  transforma¬ 
tion.  The  trailing  vortex  wake  was  modeled  by  a  sequence  of 
discrete  point  vortices,  and  the  pitching  motion  of  the  air¬ 
foil  was  modeled  by  a  series  of  small  incremental  changes 
in  angle  of  attack,  Aa,  over  a  short  period  of  time.  At. 

The  rate  of  change  of  angle  of  attack,  a,  was  then  defined 
as  Aa/At.  After  each  time  change  At,  a  was  changed  by  an 
amount  Aa .  A  discrete  vortex  was  introduced  into  the  wake  at 
a  distance  U^At  behind  the  airfoil  trailing  edge,  and  a 
bound  vortex  of  equal  strength  but  opposite  sense  was  intro¬ 
duced  to  satisfy  the  Kutta  condition  and  keep  the  total  cir¬ 
culation  in  the  flow  field  equal  to  zero.  As  each  new  vor¬ 
tex  pair  was  introduced,  all  other  trailing  vortices  were 
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assumed  to  move  in  the  wake  by  a  distance  UAt,  where  U  is 
the  velocity  induced  at  a  vortex  position  by  all  other  trail¬ 
ing  vortices,  the  bound  vortices,  and  the  free  stream  flow. 
The  unsteady  Bernoulli  equation  was  solved  using  numerical 
integration  and  differentiation  techniques  to  determine  pres¬ 
sure  difference  distribution,  vorticity  distribution,  and 
coefficient  of  lift  on  the  airfoil  for  that  instant  in  time. 
This  information  was  then  used  to  investigate  the  overall 
effect  of  constant-a  flow  as  well  as  the  effect  of  thickness 
and  camber  on  the  constant-a  problem,  and  simple  rules  for 
predicting  the  effects  were  developed. 
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THE  EFFECT  OF  TRAILING  VORTICES  ON  THE  PRODUCTION 
OF  LIFT  ON  AN  AIRFOIL  UNDERGOING  A  CONSTANT 
RATE  OF  CHANGE  OF  ANGLE  OF  ATTACK 

I .  Introduction 

It  has  been  determined  experimentally  that  an  airfoil 
pitching  at  some  rate  of  change  of  angle  of  attack  a  stalls 
at  a  higher  angle  of  attack  a  than  the  static  stall  a.  Max 
von  Kramer  first  showed  this  with  his  experiments  in  1932 

(1) ,  where  he  held  the  airfoil  fixed  in  space  and  rotated 

the  flow  over  the  airfoil  to  create  an  a.  Deekens  and  Kuebler 

(2)  and  Daley  (3)  ran  similar  experiments  for  a  constant  a, 
but  rather  than  rotating  the  flow,  they  rotated  the  airfoil 
in  a  constant  velocity  free  stream  to  produce  their  a.  In 
all  three  cases  the  stall  occurred  at  a  higher  angle  of  attack 
than  the  static-stall  angle  of  attack.  However,  because  of 
the  different  methods  used  to  produce  a,  Kramer's  results 
showed  a  much  smaller  change  in  stall  angle  of  attack  than 
did  Deekens  and  Kuebler  and  Daley. 

Following  these  experiments,  attempts  have  been  made  to 
analytically  model  the  case  of  an  airfoil  undergoing  a  con- 
stant  a.  Docken  (4)  and  Lawrence  (5)  have  tackled  the  prob¬ 
lem  using  a  momentum  integral  method,  but  both  assumed  in 
their  solution  that  the  effect  of  the  trailing  vortices  in 
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the  airfoil  wake  was  small  and  could  be  neglected.  Thus, 
they  assumed  that  the  inviscid  flow  velocity  outside  the 
airfoil  boundary  layer  at  any  angle  of  attack  was  that  which 
would  exist  in  the  steady  state  at  that  angle  of  attack.  It 
is  the  intent  of  this  thesis  to  determine  the  validity  of 
that  assumption  by  analyzing  the  effect  a  trailing  vortex 
wake  has  on  the  inviscid  flow  field  about  an  airfoil  under¬ 
going  a  constant  rate  of  change  in  angle  of  attack  (i.e., 
constant  a).  The  effect  of  the  trailing  vortex  wake  on  the 
flow  about  the  airfoil  can  be  analyzed  by  determining  how 

the  vorticity  distribution  and  pressure  difference  distribu- 

• 

tion  on  the  airfoil  develop  under  the  influence  of  the  a 
(taking  the  trailing  wake  into  account),  and  by  observing 
the  effect  of  the  a  on  the  vs.  a  curve. 
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II .  Solution  Development 


Solution  Overview 

Consider  an  airfoil  at  an  angle  of  attack  whicx:  under¬ 
goes  an  impulsively  started  motion  of  velocity  U  .  Assume 
the  airfoil  is  immersed  in  an  incompressible,  inviscid  fluid. 
Under  these  circumstances,  a  stagnation  point  of  the  flow 
would  occur  on  the  upper  surface  of  the  airfoil.  This  would 
imply  an  infinite  velocity  at  the  airfoil  trailing  edge.  It 
is  known,  however,  that  the  flow  at  the  trailing  edge  of  such 
an  airfoil  becomes  smooth  and  has  a  finite  velocity.  This 
is  known  as  the  Kutta  condition.  Imposing  the  Kutta  condi¬ 
tion  requires  the  formation  of  circulation  around  the  airfoil 
to  move  the  stagnation  point  to  the  trailing  edge.  This  cir¬ 
culation  can  be  modeled  as  a  vortex  bound  to  the  airfoil. 

The  total  circulation  in  the  flow  must  remain  equal  to  zero 
by  Kelvin's  theorem,  and  thus  circulation  in  the  opposite 
sense  is  shed  in  the  form  of  a  discrete  vortex  into  the  air¬ 
foil  wake.  The  strength  of  this  vortex  is  just  equal  and 
opposite  to  that  of  the  bound  vortex  on  the  airfoil.  The 
equal  and  opposite  strengths  of  the  bound  and  shed  vortices 
are  just  sufficient  to  satisfy  the  Kutta  condition  and 
Kelvin's  theorem. 

Thus,  when  the  airfoil  at  angle  of  attack  is  impulsively 
started,  circulation  about  the  airfoil  develops,  and  a  wake 
vortex  is  shed.  After  a  time  At,  this  shed  vortex  is  arbi¬ 
trarily  assumed  to  be  at  a  distance  U  At  from  the  trailing 
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edge  (6s21).  The  bound  vortex  and  shed  vortex  both  affect 
the  flow  about  the  airfoil,  and  their  strengths  are  such 
that  the  Kutta  condition  at  the  trailing  edge  is  satisfied. 
Knowing  the  strengths  of  these  vortices,  the  instantaneous 
values  of  circulation  about  the  airfoil,  as  well  as  the 
pressure  difference  distribution,  vorticity  distribution, 
and  coefficient  of  lift  on  the  airfoil  can  be  calculated. 
After  another  time  At,  the  shed  vortex  has  moved  further 
downstream  by  a  distance  U  At,  where  U  is  the  velocity 
at  the  shed  vortex  location  imposed  by  the  free  stream  and 
all  other  vortices,  including  the  bound  vortex.  In  cases 
other  than  the  impulsive-start  problem,  the  angle  of  attack 

may  also  have  changed  by  some  amount  equal  to  aAt,  where 

• 

a  is  the  average  rate-change  of  angle  of  attack  over  the 
given  time  period  At.  The  strength  of  the  first  shed  vortex 
remains  fixed,  and  thus  another  bound  vortex  and  shed  vortex 
must  be  introduced  to  keep  the  Kutta  condition  satisfied. 

This  second  shed  vortex  is  assumed  to  be  at  a  distance 
U^At  behind  the  trailing  edge.  The  equal  and  opposite 
strengths  of  these  new  bound  and  shed  vortices  are  again  de¬ 
termined  by  imposing  the  Kutta  condition.  Now,  for  this  new 
instant  in  time,  the  instantaneous  values  of  airfoil  cir¬ 
culation,  pressure  difference  distribution,  vorticity  dis¬ 
tribution,  and  coefficient  of  lift  can  once  again  be  cal¬ 
culated.  This  process  can  be  repeated  for  any  number  of 
discrete  time  steps  At  desired,  and  for  that  matter  any 
a(t),  although  in  this  study  a  was  held  constant.  By 


following  this  method,  a  time  history  of  the  development  of 
circulation,  pressure  difference  distribution,  vorticity 
distribution,  and  coefficient  of  lift  on  the  airfoil  can  be 
observed . 

Equations  for  Flow  About  a  Cylinder 

When  solving  a  problem  in  two-dimensional  incompressi¬ 
ble,  irrotational  flow,  it  is  often  useful  to  make  use  of 
conformal  mapping.  In  this  case,  the  problem  is  solved  for 
flow  about  a  two-dimensional  cylinder,  then  the  Joukowski 
transformation  is  used  to  find  the  solution  for  a  Joukowski 
airfoil . 

Consider  the  flow  of  an  incompressible,  irrotational 
fluid  in  the  o-plane.  The  flow  is  inclined  at  an  angle  a 
to  the  x-axis  (see  Fig.  1).  The  stream  function  t//  and  poten 
tial  function  0  for  this  flow  are  given  by  (7:245): 


=  U^tY  cos  a  -  x  sin  a)  (1) 

0  =  U  (X  cos  a  +  Y  sin  a)  (2) 

oo 


where  U  is  the  magnitude  of  the  free  stream  velocity. 

If  a  doublet  of  strength  K,  axis  inclined  at  angle  a 
to  the  X-axis,  is  placed  at  the  origin  of  the  o-plane,  the 
stream  function  \}/  and  potential  function  0  are  given  by: 


Y  cos  a  -  X  sin  a 

2  2 
X  +  Y 


(3) 
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X  cos  g  +  Y  sin  a 

2  2 
X  +  Y 


(4) 


Since  stream  functions  and  potential  functions  are  lin¬ 
ear,  they  may  be  superimposed  to  create  new  flows.  Therefore, 
the  stream  and  potential  functions  for'  a  doublet  in  a  uniform 
free  stream  at  angle  a  to  the  X-axis  can  be  written  as: 

U/  =  U  (Y  cos  a  -  X  sin  a)  -  ^ 

oo  Z'TT 


Y  cos  a  -  X  sin  a 

2  2 
X*  +  Y 


=  (Y  cos  a  -  X  sin  a) 


K 

2rr 


(5) 


0  =  U  (X  cos  a  +  Y  sin  a) 

OO 


X  cos  a  +  Y  sin  a 
2  2 


=  (X  cos  a  +  Y  sin  a) 


u 


_K 

2it 


X2  + 


(6) 


.  2 

If  strength  K  is  such  that  K/ZrrU^  =  a  ,  where  a  is  the 
radius  of  a  cylinder,  then  the  zero  streamline  will  be  the 
line  along  which  Y  cos  a  =  X  sin  a  and  the  cylinder  of 
radius  a  ,  centered  at  the  origin  (8:89).  The  stream  func¬ 
tion  and  potential  function  are  now: 


if/  =  U  (Y  cos  a  -  X  sin  a) 

OO 


1 


(7) 


0  =  U^X  cos  a  +  Y  sin  a) 


(8) 
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After  using  the  trigonometric  identity  (8:86) 


manipulation  of  Eq.  (18).  The  velocities  in  the  X-direction 


(Up)  and  the  Y-direction  (V^)  ares 


0  '  3Y 


U 


cos  a 


2  2  2  2 
a  cos  a  _  2a  XY  sin  a  2a  Y  cos  a 

x2  +  Y2  (x2  +  Y2)2  (x2  +  Y2)2 


N 

+  E 
i=l 


a 

2tt 


Y-B. 

_ l _ 

(X-Ai)2+(Y-Bi)2 


Y-rJi 


(X-  ^i)2+Y-TJi)2 


(20) 


V 

0 


-a^ 

ax 


u 


sin  a 


sin  a 


2  2 
X  +  Y 


2  2  2 
2a  XY  cos  a  2a  X  sin  a 

(X2  +  Y2)2  (X2  +  Y2)2 


N 
+  E 
i=  1 


a 

2rr 


X-^i 

(X-^)2+(Y-tj.)2 


X-A. 

_ l _ 

(X-Ai)2+(Y-Bi)2 


(21) 


The  stream  function  ,  potential  function  0,  and  general 
velocities  UQ  and  for  any  point  in  the  0-plane  are 

now  known . 


Joukowski  Transformation 

The  Joukowski  transformation  can  be  used  to  transform  a 

flow  from  a  cylinder  plane  to  an  airfoil  plane.  We  already 

have  the  equations  for  the  flow  in  the  0-plane,  where  the 

cylinder  is  centered  at  the  origin.  In  complex  variable 

form,  the  position  of  a  point  in  the  0-plane  can  be  expressed 
i  0 

as  0  =  r  e  .  We  can  transform  this  point  to  the  0’ -plane 

—  i8 

by  the  transformation  0 '  =  oe  p  +  (i  (10:461).  This  ro¬ 


tates  points  in  the  0-plane  8  radians  clockwise,  and  then 


displaces  them  by  an  amount  4.  The  positive  X-coordinate 

crossed  by  the  cylinder  in  the  c-plane,  p  ,  maps  into  the 

point  Ct '  (see  Fig.  2).  Transforming  from  the  p’-plane 

to  the  Z-plane,  the  Joukowski  transformation  is  used.  It 

is  given  by:  ,2 

Z  =  0*  +  ~T~  (22) 

This  transforms  the  cylinder  in  the  o'-plane  to  an  airfoil 
shape  in  the  Z-plane.  The  point  p  •  maps  to  the  trailing 
edge  of  the  airfoil  shape  in  the  Z-plane. 

Determination  of  S trenqth  of  Vortices 

Having  seen  how  the  Joukowski  transformation  maps  a 
cylinder  in  the  o-plane  into  an  airfoil  shape  in  the  Z-plane, 
it  is  time  to  relate  the  flows  in  the  0  and  Z-planes.  As 
mentioned  before,  the  solution  involves  placing  discrete 
vortices  in  the  airfoil  wake  to  simulate  the  vortex  sheet 
shed  into  the  wake  as  circulation  builds  around  the  airfoil. 
Images  of  these  vortices  are  placed  inside  the  cylinder  to 
simulate  the  airfoil  bound  vortex.  The  strength  of  each 
vortex  is  determined  by  satisfying  the  Kutta  condition  at 
each  discrete  time  step.  The  Kutta  condition  implies  that  a 
stagnation  point  of  the  flow  is  at  the  airfoil  trailing  edge. 
As  seen  in  the  Joukowski  transformation,  the  point  ot  in 
the  p-plane  maps  into  the  trailing  edge  of  the  airfoil  in  the 
Z-plane.  Therefore,  establishing  a  stagnation  point  at  p 
in  the  o-plane  satisfies  the  Kutta  condition  in  the  Z-plane 
(10:469) . 
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p '  -  Plane 
Y' 


of  the  Joukowski  Transformation 
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To  make  ofc  a  stagnation  point  in  the  o-plane,  set  the 
velocity  Un  equal  to  zero  at  that  point  and  solve  for  the 
circulation  strength  Pin  Eq.  (20).  This  value  of  P  will 
be  the  same  in  the  2-plane,  for  circulation  is  unchanged  in 
the  Joukowski  transformation  (10j458).  Taking  Eq.  (20)  for 

Up,  letting  the  cylinder  radius  a  =  1  ,  and  recalling  that 

2  2  .  . 

X  +  Y  =1  for  points  on  a  cylinder  of  radius  a  =  1  ,  one 

arrives  at: 


U  =  0  =  U  (cos  a  -  cos  a  -  2XY  sin  a  +  2Y2  cos  a 
o  00 1 


N  Pi 

+  I  -=-=■ 

i-1  21T 


Y-B. 

l 


Y-\ 


L(X-Ai)2+(Y-Bi)2  (X-^)2+(Y-l?i)2J 


(23) 


Let  N  =  1  to  solve  for  the  strength  of  the  first  shed  vor¬ 
tex.  Solving  for  P/U .  one  gets: 


2tt  ( 2XY  sin  a  -  2Y2  cos  a) 

Y-B _ Y-  V 

(X-A)2+(Y-B)2  (X-^)2+(Y-n)2 


(24) 


Recalling  the  relations  for  A  and  B  given  by  Eqs .  (13)  and 

2  2 

(14)  and  again  making  use  of  the  fact  that  X  +  Y  -  1  on 
the  cylinder,  Eq.  (24)  becomes 


r_  a  4tt  (X  sin  a  -  Y  cos  a)(52  +  rj2  +  i  -  2(xg+Y»7)) 


u 


»2  - 


(25) 


★  , 

Define  r/ U  as  P  ,  the  non-dimensional  circulation,  and 

00 

note  that  at  ot»  X  =  1  and  Y  =  0  .  This  changes  Eq .  (25) 
to 
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(26) 


Velocities  Induced  at  Discrete  Vortices  and  on  the  Cylinder 
For  each  time  step  taken,  the  strength  of  the  vortex 
pair  introduced  at  that  time  step  can  now  be  calculated 
using  Eq.  (27).  However,  from  one  time  step  to  the  next, 
each  vortex  introduced  in  the  wake  moves  away  from  the  air¬ 
foil  some  distance,  that  distance  being  equal  to  the  velocity 
at  the  position  of  the  vortex  times  the  time  step.  At  . 

The  velocity  at  the  position  of  each  vortex  depends  upon  not 
only  the  free  stream  velocity,  but  also  the  velocities  in¬ 
duced  at  that  position  by  all  other  vortices  in  the  field. 
Equations  (20)  and  (21)  can  be  used  to  find  that  velocity  in 
the  o-plane.  Let  (t^^L,)  be  the  coordinates  of  the  posi- 
tion  at  which  the  trailing  vortex  is  located,  the  i  sub¬ 
script  denote  the  time  step  at  which  the  velocity  is  com¬ 
puted,  and  the  k  subscript  identify  an  individual  vortex. 
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Solving  equations  (20)  and  (21)  for  the  non-dimensional 

U  /U  ,  V  /U  ,  and  recalling  P  =  P  /U  ,  one  gets; 

m  ■»  m  30  a  m  m  °°  * 


U  .  2 

m  cos  a.  -  cos  a.  -  2£  t)  sm  a.  +  2n  cos  a- 
rr-  =  l  _ l  ^m  m _ i  m _ i__ 

00  .  2  _  2  ,  A  2  ^  2N  2  ,  A  2  ^  2  N  2 


I  +  T1 

’m  m 


U  "  +  n  ")  ) 

m  m  mm 


r  * 

*m 

2tt 


»JmU  2  +  n_2  ~  1) 

'  m _ m _ m 


|(  £  2  +17  2)2  -  2  £  2  -  2  T)  2  +  lj 

L  x  m  m  J 


(28) 


k  =  l, 
k/m 


Ik 

2tt 


^m_Bk 


(A  -A.  )2+(r)  -B,  )2  (£  )2+(tj  -ti  )2 

sm  k  ’m  k  sn  sk  'm  ’k 


Vm  _  sin  a-  -  sin  a.  -  2£  tj  cos  ct.  +  2£  2  sin  a- 
—  1  1  mm  1  m  1 


U 


$  2+T)  2  (£  2  +  n  2)2  U  2  +  rj  2)2 

sm  'm  vm  m  m  ' m 


'm 


2tt 


*  "m2  ~  11 


($  2  +  n  2)  -  2£  2  -  2?7  2  +  1 

sm  m  sm  'm 


(29) 


N 


rK 


I  2tt 
k=l, 
k/m 


^m  ^k 


t(<m-«K)2+<V\>2  <«m-AK,2t("n-BK>2 
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Equations  (28)  and  (29)  describe  the  velocity  that  exists  at 


any  vortex  location  ( £  ,  V  )  as  induced  by  the  uniform 

free  stream  and  all  other  vortices  in  the  flow  field. 

The  velocity  induced  along  the  cylinder  surface  is 

needed  to  determine  the  pressure  distribution  on  the  airfoil. 

Again  making  use  of  the  o-plane  velocity  equations  (20)  and 

2  2 

(21),  and  recalling  that  X  +  Y  =1  on  the  cylinder  in 
the  0-plane,  the  non-dimensional  velocities  become: 


=  -2XY  sin  a 


2Y‘ 


N  j-.* 

cos  a  +  I  *i 
i  =  l  2-tt 


Y(^i2+^i2-l) 


UiH2-2(xVYV+1 


(30) 


N  * 

V  _  -2XY  cos  a  +  2X2  sin  a  -  E  ^i 
U  “  i  =  l  2rr 


X(^i2+7?i2_1) 


h2*\ 2-2<ni*T'vi)+i 


(31) 

where  (X,Y)  are  coordinates  of  a  point  on  the  cylinder  sur¬ 
face.  In  the  c-plane,  these  points  can  be  easily  put  into 
cylindrical  coordinates.  Let  X  =  r  cos  9,  Y  =  r  sin  9  , 

where  r  =  1  ,  the  radius  of  the  cylinder.  Equations  (30) 

and  (31)  then  become: 

U  2 

—  =  -2  cos  0s  in  6  sin  ci  +  2  sin  9  cos  a 

OO 

sin  9  (  ^2  +  ?7i2  -  1 ) 

^2  +  V ^2  -  2(£^  cos  0  +  ^  sin  9)  +  1  . 

(32) 


r 

+  E  'i 
i  =  1  2tt 
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—  =  -2  cos  0  sin  0  cos  a  +  2  cos  6  sin  a 


?  r* 

+  E  i 
i  =  l  2tt 


cos  0  ( I i2  +  *7j_2  -  1) 


L  t,2  +  -  2 ( ^ i  cos  0  +  /7i  sin  0) 


1  J 

(33) 


Since  the  cylinder  surface  is  a  streamline  of  the  flow, 
the  velocity  on  the  cylinder  is  always  parallel  to  the  sur¬ 
face.  Therefore,  the  magnitude  of  the  non-dimensional  ve¬ 
locity  tangent  to  the  surface  of  the  cylinder,  U6  ,  is 
justs 


U0  = 


(34) 


Circulation  About  the  Airfoil 

As  the  wake  behind  the  airfoil  forms,  circulation  de¬ 
velops  about  the  airfoil  in  the  form  of  a  bound  vortex. 

The  strength  of  this  bound  vortex  defines  the  total  circula¬ 
tion  about  the  airfoil.  Since  the  value  of  the  total  cir¬ 
culation  in  the  flow  field  must  be  zero,  then  the  strength 
of  the  circulation  about  the  airfoil  must  be  equal  in  magni¬ 
tude  and  opposite  in  sign  to  the  total  circulation  in  the 
wake.  The  circulation  in  the  wake  is  just  the  sum  of  the 
strengths  of  all  the  discrete  vortices  in  the  wake.  The  cir¬ 
culation  can  be  calculated  by 
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where  r  is  the  non-dimensional  circulation  about  the  air- 

cl 

*  .  . 

foil,  77  is  as  defined  in  Eq .  (27),  and  N  is  the  number  of 
discrete  vortices  in  the  wake. 

Velocity  in  the  Airfoil  Frame 

Now  that  expressions  for  velocities  in  the  o-plane  are 
known  (Eqs.  (32),  (33)  and  (34)),  their  values  at  correspond¬ 
ing  points  in  the  Z -plane  can  be  found.  Consider  the  complex 
potential  F(o)  =  0  +  i ty,  where  0  =  X  +  iY  in  the  o-plane. 

The  complex  velocity  in  the  o-plane  is  dF/do  =  w(o)  . 

—  i  B 

Since  0 '  =  oe  +  4  ,  then  0 '  is  a  function  of  0 •  By 
the  chain  rule  of  differentiation,  dF/do  =  (dF/dc * ) (dp ’/do)  , 
where  dF/do'  =  w(o')  »  the  complex  velocity  in  the  o' -plane. 
Therefore 

w(O')  =  w(o)  .  _L_  .  (36) 

do  ' 
do 

Note  that  the  magnitude  of  do’/do  =  1  ,  and  thus  the  magni¬ 
tude  of  the  complex  velocity  in  the  o-plane  equals  the  mag¬ 
nitude  of  the  complex  velocity  at  the  corresponding  point  in 
the  o' -plane.  By  similar  arguments,  knowing  that  the  trans¬ 
formation  from  the  o' -plane  to  the  Z-plane  is  given  by  Eq. 
(22),  and  using  the  chain  rule  once  more,  it  can  be  shown 
that 

U„  -  iV_  =  w(Z)  =  w(o')  .  _L_  =  w(o)  .  _±_  .  (37) 

Z  Z  dZ_  dZ 

do ’  do ' 
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Thus,  given  the  magnitude  of  a  velocity  at  a  point  in  the 

O-plane,  the  magnitude  of  the  velocity  at  the  corresponding 

point  in  the  Z-plane  can  be  found  using  Eq .  (37). 

All  the  tools  needed  to  analyze  the  airfoil  wake  effects 

are  now  known.  The  stream  and  potential  functions  in  the 

O-plane  are  known,  from  which  velocities  in  the  o-plane  can 

_  * 

be  found.  The  strength  1  of  the  wake  vortices  and  their 

images  can  be  found  by  requiring  the  stagnation  point  in  the 

O-plane  to  remain  fixed,  thus  satisfying  the  Kutta  condition 

_* 

on  the  airfoil  in  the  Z-plane.  The  values  of  l  are  the 
same  in  both  planes,  and  velocities  in  the  o-plane  can  be 
directly  transformed  to  the  Z-plane. 

Pressure.  Lif t .  Vorticity  Distribution  on  the  Airfoil 

The  unsteady  Bernoulli  equation  is  used  to  calculate 
the  pressure  on  the  airfoil.  It  is  given  by  (6:18) 

P+i5  0U2  +  o|f=D  f(t)  (39) 

where  P  is  pressure,  o  is  the  fluid  density,  U  is  the  fluid 
velocity,  0  is  the  potential  function,  and  f(t)  is  a  func¬ 
tion  of  time  independent  of  position.  The  subscripts  i  and 
u  will  be  used  to  denote  the  lower  and  upper  surfaces  of  the 
airfoil,  respectively.  Subtracting  the  pressure  on  the  upper 
surface  from  the  pressure  on  the  lower  surface  yields 


(40) 


t 


J' 


The  velocities  in  Eq .  (40)  are  those  tangent  to  the  airfoil 
surface.  Along  a  streamline,  U  =  30/3S  ,  where  S  is  the 

coordinate  along  the  streamline.  This  implies  that 
0  =  /  UdS  .  Therefore, 


p,  -  P, ,  -  %  0<u„2  -  u  2)  ♦  ofj-  fSJ  (Ct  -  u.)  dS 


u 


U 


3 1 


'u 


(41) 


Integrating  in  Eq.  (41)  from  to  S  ,  the  point  on  the 

a 


U  and  U.  are  known  or  desired, 

u  Z 


streamline  where  P 
one  finds  that  the  integral  is  zero  from  to  -He  , 

the  airfoil  leading  edge.  Equation  (41)  then  becomes 


3  r'Sa 


pe  -  pu  =  *  0(uu2  -  u/>  +  Uc(uu  -  V 


(42) 


Introduce  the  following  non-dimensional  variables,  identified 
by  the  superscript  *: 


*  t  U 

★  oo 


t  =  r 


>  u  =  §-  ■  s  “  fe  <43> 

OO 


Equation  (42)  now  becomes: 


*  ★ 


U  *  pu  *  ^.2(,JU*2  -  u/2)  *  °u„2i_  /.*  (Uu  -up)  ds  . 


3 1 


(44) 


Recall  the  definition  of  the  coefficient  of  pressure 


C  =  P-P  /H  CU  .  The  difference  between  the  coefficients 


OO  oo 


of  pressure  on  the  lower  and  upper  airfoil  surface  is  called 


ACp  .  Using  this  definition,  one  obtains  from  Eq.  (44) 
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P£-Pu  3  „Sa  *  *  * 

-JL— f  =  -  AC  =  U  -  U  +  2  —  S_l  (U  -U  )  dS  .  (45) 

i  z  pu  x  at  1  u  1 

oo 

The  coefficient  of  lift  per  unit  span,  ,  is  de- 

fined  as  =  L/hpU ^  c  ,  where  L  is  the  lift  per  unit  span 

and  c  is  the  chord.  Since  L  is  defined  as  f  fcAPdX  ,  and 

-He 

AP  =  -ACp  •  ho^J'  ,  then 

r  -  ACP  dx  (46) 

Ci  "  J -^c  c 

where  X  is  measured  along  the  chord  of  the  airfoil. 

★ 

The  non-dimensional  vorticity,  y  ,  is  defined  as 
★  ★ 

Uu  -  .  The  vorticity  distribution  can  be  easily  cal- 

★ 

culated  using  this  relation,  since  the  value  of  y  can  be 
obtained  directly  for  any  position  along  the  airfoil  chord 
where  the  velocities  on  the  upper  and  lower  surfaces  of  the 
airfoil  are  known. 


Numerical  Solution  Process 

The  following  procedure  is  used  to  numerically  analyze 
the  wake  vortex  effects  on  the  airfoil. 

★ 

Step  1  -  Select  a  non-dimensional  time  step  At  ,  defined 

as  At  =  AtU^/ ( Hc)  •  Let  the  airfoil  begin  its  motion  ,it 

an  initial  angle  of  attack  aQ  and  velocity  .  Assume 

★ 

that  after  a  time  At  ,  at  time  i  -  1  ,  the  first  shed 
vortex  is  at  a  position  U At  downstream  of  the  airfoil 


trailing  edge  in  the  direction  of  the  velocity  at  the  trail¬ 


ing  edge. 


Step  2  -  The  first  vortex  is  thus  at  (|,77)  in  the  o-plane, 

* 

and  Eq .  (26)  can  be  solved  for  /  .  The  circulation  about 

.  •  * 
the  airfoil  is  -p  . 

Step  3  -  Eqs .  (32),  (33)  and  (34)  can  now  be  solved  for  the 
velocity  in  the  p-plane  at  any  point  (X,Y)  on  the  cylinder. 
These  velocities  can  be  transformed  to  the  Z-plane  by  Eq.  (37) 
Step  4  -  Eq.  (45)  can  be  solved  for  AC  using  a  trapezoidal 

r"' 

rule  with  variable  AX  for  the  integration  along  the  upper  and 
lower  surfaces  of  the  airfoil,  and  a  three-point  backward 
difference  differentiation  approximation  for  the  derivative 
with  respect  to  time.  (For  steps  i  =  1  and  i  =  2  ,  a  two- 
point  linear  difference  method  is  used  for  the  time  deriva¬ 
tive.)  Eq .  (46)  can  then  be  solved  for  ,  again  using  a 
trapezoidal  rule  with  variable  AX  for  the  integration.  The 
non-dimensional  vorticity  distribution  can  be  calculated 

*  *  if 

directly  as  y  =  uu  ~ 

Step  5  -  The  velocity  of  the  shed  vortex  is  calculated  using 

Eqs.  (28),  (29)  and  (37).  For  the  next  time  step,  the  vor- 

★  ★ 

tex  has  moved  in  the  Z-plane  by  an  amount  Usv  At  ,  where 

★ 

Usv  is  the  non-dimensional  velocity  of  the  vortex  just  cal¬ 
culated.  Its  position  in  the  o-plane  is  then  determined  by 
the  inverse  Joukowski  transformation,  given  by 

Z  ±  (Z2  -  4ofc2 )  ^  -  ii 
2 

where  only  the  plus  sign  of  the  ±  term  gives  a  value  of  o  in 
the  wake,  and  thus  it  is  the  value  used  for  c. 
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i  =  2  ,  the  vortex  shed 


at  time  i  =  1  is  at  the  position  computed  in  step  5.  As¬ 
sume  the  vortex  shed  at  time  i  =  2  is  at  the  position 
U^At  downstream  of  the  trailing  edge.  The  angle  of  attack 

at  l  =  2  is  now  a  +  a  At  ,  where  a  is  defined  as 

o 


•  * 

a 


a 


u 


(48) 


Step  7  -  All  terms  in  Eq.  (27)  are  now  known,  and  this  equa- 

★ 

tion  can  be  solved  for  /T,  .  The  circulation  about  the  air¬ 
foil  is  -(  /J*  +  )  • 

Step  8  -  Eqs .  (32),  (33)  and  (34)  can  be  solved  for  the  ve¬ 
locity  on  the  cylinder  in  the  o-plane,  and  then  the  veloci¬ 
ties  can  be  expressed  in  the  Z-plane  using  Eq.  (37). 

Step  9  -  The  values  of  AC^  and  can  be  found  using  Eqs. 

(45)  and  (46)  respectively,  and  the  vorticity  distribution 

*  *  -fc 

Y  is  again  just  Uu  -  . 

Step  10  -  The  velocity  at  the  shed  vortices  can  be  calculated 

using  Eqs.  (28),  (29)  and  (37).  The  vortices  are  then  moved 

★  ★  , 

in  the  Z-plane  a  distance  Usv  At  .  This  new  position  of 
each  vortex  in  the  Z-plane  is  the  assumed  location  of  each 
vortex  for  the  next  time  step.  Each  vortex  position  in  the 
0-plane  can  be  determined  by  using  the  inverse  Joukowski 
transformation,  Eq .  (47). 

Step  11  -  For  each  time  step  i,  the  position  of  each  vortex 
is  known  from  time  step  i  -  1,  and  the  vortex  shed  at  time 
i  is  assumed  to  be  a  distance  U  At  downstream  of  the  trail- 

OO 


ing  edge  in  the  direction  of  the  velocity  at  the 


trailing  edge.  The  angle  of  attack  a  at  time  step  i  is 
ai  =  ai-i  +  <**At*  •  Eqs-  (27),  (35),  (32),  (33),  (34),  (37) 

(45)  and  (46)  are  then  used  to  compute  [ 7*  ,  airfoil  cir¬ 

culation,  velocities  on  the  airfoil,  AC  and  C»  for  time 

P  * 

step  i.  Eqs .  (28),  (29)  and  (37)  are  used  to  find  the  posi¬ 
tion  of  the  shed  vortices  for  time  step  i  +  1. 

Step  11  is  repeated  as  often  as  desired  to  compute  the 
airfoil  circulation,  pressure  difference  distribution,  coef¬ 
ficient  of  lift,  y  distribution,  and  shed  vortex  positions 
for  any  discrete  time  desired. 


Ill .  Results 


Numerical  Method  Ver if ication 

Before  exploring  the  effect  of  a  constant-a  flow  on  the 
production  of  lift  on  a  Joukowski  airfoil,  the  method  devel¬ 
oped  here  was  compared  to  the  results  of  others.  The  first 
test  case  was  that  of  a  flat  plate  impulsively  started  at  an 
infinitesimal  angle  of  attack,  a.  This  problem  was  first 
explored  by  Wagner  (11)  in  1925.  Wagner  assumed  in  his  analy¬ 
sis  that  the  wake  vortex  sheet  remained  along  the  x  -axis 

d. 

at  all  times.  For  infinitesimal  a,  this  is  a  good  approxi¬ 
mation.  In  Figs.  3  and  4  a  numerical  computation  for  a  flat 

plate  in  the  Z-plane,  impulsively  started  at  a  =  0.01 
★ 

radians.  At  =0.02  ,  is  compared  with  Wagner’s  analytic 

results  and  Giesing's  (12)  numerical  results.  The  horizon¬ 
tal  axis  scale  is  U  At/^c,  which  is  the  non-dimensional  dis- 
tance  the  airfoil  has  traveled  since  the  motion  started, 
having  a  value  of  one  for  each  half-chord  distance  of  air¬ 
foil  translation.  The  vertical  scales,  r/rss  and  Cz/Cz  , 

ss 

are  the  ratios  of  F  or  C ^  to  the  steady-state  values  that 
would  be  obtained  after  a  long  period  of  time  has  passed, 
respectively.  Both  the  build-up  of  circulation,  r  (Fig.  3), 
and  the  coefficient  of  lift,  (Fig.  4),  closely  approximate 
Wagner's  curves . 

In  1977,  Shung  (6)  developed  a  numerical  method  similar 
to  the  one  presented  in  this  thesis,  but  limited  the  study 

*■ 

*  to  that  for  a  flat  plate  at  a  constant  a.  Unlike  Wagner, 
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however,  Shung's  method  allowed  for  vortex  interactions  in 
the  airfoil  wake,  as  does  the  method  in  this  thesis.  The 
major  difference  between  Shung's  numerical  technique  and  the 
technique  presented  here  lies  in  the  method  used  to  integrate 
the  Unsteady  Bernoulli  Equation,  Eq .  (44).  Since  Shung  was 
limited  to  a  flat  plate,  he  was  able  to  integrate  using  the 
Gauss-Chebyshev  quadrature  formula  (6:21).  In  this  thesis, 
a  trapezoidal  rule  with  variable  AX  was  used  for  that  inte¬ 
gration.  This  allowed  easy  application  to  Joukowski  air¬ 
foils.  Figure  5  depicts  C„/C^  for  a  flat  plate  at 

ss  * 

a  =0.1  radians  at  various  values  of  At  .  Comparing  the 

numerical  solution  with  Wagner's  curve,  one  sees  that  for 

★ 

smaller  values  of  At  ,  the  numerical  solution  approaches 
Wagner's  analytic  solution.  Shung  (6:45)  noted  the  same 

tendency  with  his  numerical  solution.  Figure  6  depicts  the 

★ 

formation  of  the  vorticity  distribution  y  on  a  flat  plate 

★ 

at  a  =  0.1  radians.  At  =0.02  .  Note  that  immediately 

after  the  airfoil  begins  its  motion,  the  vorticity  is  nega¬ 
tive  near  the  trailing  edge,  but  as  time  passes,  the  vortic¬ 
ity  distribution  approaches  that  for  the  steady-state  condi¬ 
tion.  Shung  (6:44)  showed  the  same  effect  in  his  study 

★ 

using  a  =  0.1  radians.  At  =  0.1  .  In  fact,  his  results 

are  identical  to  the  results  presented  in  Fig.  6.  Similarly, 
the  build-up  to  the  steady-state  pressure  difference  distribu¬ 
tion  for  a  flat  plate  at  a  =  0.1  radians  can  be  seen  in 
Fig.  7.  Shung  also  demonstrated  wake  vortex  sheet  roll-up 
behind  an  impulsively-started  flat  plate.  The  method  of 
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this  thesis  demonstrates  the  same  phenomenon,  as  can  be 
seen  in  Fig.  8,  which  depicts  wake  vortex  sheet  roll-up  for 

the  case  of  an  impulsively-started  flat  plate  at  a  =  0.1 

★ 

radians,  At  =0.1  ,  the  same  conditions  depicted  by 

Shung  (6:43).  Shung's  depiction  and  Fig.  8  are  almost  iden¬ 
tical  . 

The  flat  plate  is,  in  fact,  a  special  case  of  the  more 

general  Joukowski  airfoil,  and  other  airfoils  in  this  family 

have  been  studied.  Giesing  (12)  also  developed  a  numerical 

procedure  to  account  for  wake  effects  on  the  build-up  of 

lift  on  an  arbitrary  airfoil.  Giesing  published  a  curve  of 

C./C  for  a  25.5%  thick  symmetric  Joukowski  airfoil  impul- 
1  ss 

sively  started  at  a  =  0.01  radians.  Using  the  same  airfoil 

and  motion  conditions,  a  C ^/C ^  curve  was  developed  using 

ss 

the  numerical  method  presented  in  this  thesis.  Figure  9 
compares  those  two  curves  with  Wagner's  curve  for  a  flat 
plate.  As  can  be  seen  in  Fig.  4,  Giesing's  curve  predicts 
values  below  those  predicted  by  the  present  method.  Whereas 
the  present  method  over-predicted  for  a  flat  plate,  and 
Giesing's  method  under-predicted  for  a  flat  plate,  it  seems 
likely  that  the  ideal  solution  is  bracketed  by  the  present 
method  and  Giesing's  method.  Both  curves  for  the  25.5% 
thick  symmetric  Joukowski  airfoil  show  a  greater  delay  in 
lift  production  than  does  the  flat-plate-airfoil  curve.  This 
agrees  with  an  analysis  done  by  Chow,  who  showed  that  air¬ 
foils  of  increased  thickness  develop  lift  at  a  slower  rate 
than  thinner  airfoils  (13:14). 


Vortex  Sheet  Roll-up  For  a  Flat  Plate,  a=0.1  Radians 


Thick  Symmetric  Joukowski  Airfoil 


Giesing's  numerical  technique  and  the  one  presented  in 
this  thesis  differ  in  several  ways.  One  major  difference 
between  the  two  numerical  techniques  is  in  the  way  the  motion 
of  the  discrete  vortices  in  the  airfoil  wake  is  predicted. 

In  Giesing's  technique,  after  each  time  step,  the  non-dimen¬ 
sional  velocity  induced  at  each  trailing  vortex  position, 

*  * 

,  is  calculated  and  then  multiplied  by  At  to  approximately 

predict  where  that  discrete  vortex  will  be  at  the  next  time 

step.  The  non-dimensional  velocity  induced  at  that  pre- 
«  .  ★ 

dieted  position,  (Jc  ,  is  then  calculated.  The  average  of 

*  ★  ,  ★ 

+  Uc  is  then  multiplied  by  At  to  correct  the  predic¬ 
ted  position  of  each  discrete  vortex  for  the  next  time  step. 
This  method  can  be  referred  to  as  a  Predictor-Corrector 
method.  The  numerical  technique  presented  in  this  thesis  pre¬ 
dicts  the  discrete  vortex  position  in  the  same  manner  as 
Giesing's  predictor,  but  no  corrector  velocity  is  computed 
or  used.  The  predicted  velocity  is  the  only  velocity  used 
to  update  vortex  position. 

To  determine  the  effect  a  Predictor-Corrector  method 
has  on  the  numerical  solution,  a  program  incorporating 
Giesing's  Predictor-Corrector  method  was  developed.  The  re¬ 
sults  obtained  using  this  program  were  compared  with  the  re¬ 
sults  obtained  using  the  method  presented  in  this  thesis 
for  the  same  airfoil  and  conditions  of  motion.  Figure  10 
shows  a  comparison  of  wake  shape  as  computed  by  the  two 
methods  for  an  impulsively  started  flat  plate  airfoil  at 

O 

a  =  10  .  Only  in  the  area  of  starting  vortex  roll-up 


does  the  difference  in  vortex  position  become  apparent. 
Further,  the  values  of  C^/C^  computed  using  the  two  methods 
are  nearly  identical,  as  seen  in  Table  I.  A  comparison  be¬ 
tween  the  two  methods  was  also  made  for  the  case  of  the  flat 

•  •  0  •  ★ 
plate  initially  at  a  =  0  subjected  to  an  a  =  0.035 

One  can  see  in  Table  II  that  the  values  of  C^/C^  computed 

ss 

using  the  two  methods  are  once  again  nearly  identical.  It 

was  thus  determined  that  the  added  computation  time  incurred 

by  using  the  Predictor-Corrector  method  was  not  needed,  and 

thus  not  included  in  other  studies  in  this  thesis. 

As  a  final  check  of  the  method  of  this  thesis,  the  de- 
★ 

velopment  of  y  and  AC^  on  an  impulsively-started  symmetric 
airfoil  with  thickness  was  determined.  Figures  11  and  12 
show  that  for  a  25.5%  thick  symmetric  Joukowski  airfoil  im- 

•k  % 

pulsively  started  at  a  =  0.1  radians,  y  and  AC  build  to 

P 

their  steady-state  values  in  much  the  same  manner  as  for  a 
flat  plate.  Figs.  6  and  7,  when  using  the  numerical  method 
of  this  thesis. 

It  has  been  shown  that  for  the  test  cases  above,  the 
numerical  method  presented  here  is  in  good  agreement  with 
the  work  of  others  (6  ;  1 1 j 1 2 ; 13 )  .  There  are  two  major  advan¬ 
tages  in  using  this  numerical  method  over  other  methods. 
First,  unlike  Shung,  one  is  not  limited  to  a  flat  plate. 
Second,  comparing  to  Giesing,  the  simpler  method  of  vortex 
motion  prediction  greatly  decreases  computer  run  time  while 
having  a  negligible  effect  on  the  prediction  of  lift  build¬ 
up  on  an  airfoil. 
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TABLE  I 


Comparison  of  C,  Calculated  by  Simple-Predictor 
Method  to  C ,  Calculated  by  Predictor  Corrector 
Method.*  Flat  Plate  Airfoil,  a  =  10 


t* 

C2 

Predict or -Cor rector 

C* 

1 

.70677 

.70779 

2 

.75962 

.75974 

3 

.79825 

.79827 

4 

.83271 

.83272 

5 

.85997 

.86000 

6 

.88199 

.88203 

7 

.90008 

.90010 

8 

.91515 

.91518 

9 

.92788 

.92789 

10 

.93873 

.93875 

TABLE  II 

Comparison  of  C^  Calculated  by  S imple- Predictor 
Method  to  C.  Calculated  by  Predictor  Corrector 
Method.  Flat  Plate  Airfoil,  a*  =  0.035 


★ 

t 

CZ 

Pr edict or-Cor rector 

0.2 

.  13710 

.13968 

0.4 

.16284 

.16543 

0.6 

.18888 

.19152 

o 

• 

00 

.21544 

.21814 

1 .0 

.23061 

.23062 

1.2 

.27020 

.27301 

1.4 

.29837 

.30122 

1.6 

.32702 

.32991 

1.8 

.35612 

.35906 

2.0 

.38564 

.38863 
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Application  to  Constant-a  Flow 

In  the  last  section  it  was  shown  that  the  present  method 
compares  favorably  to  the  results  of  others  for  the  con¬ 
stant-a,  impulsively-started  airfoil.  In  this  section  the 
results  of  applying  the  method  to  the  previousiy-unstudied 
problem  of  constant-a  flow  is  presented.  The  presentation 
of  these  results  is  broken  into  four  parts  in  order  to  more 
systematically  explore  and  understand  the  interplay  of  pos¬ 
sible  effects.  These  four  parts  deal  with  the  effects  and 
selection  of  starting  conditions,  the  general  effect  of  a 
on  the  build-up  of  C^,  the  effect  of  thickness,  and  the  ef¬ 
fect  of  camber,  respectively. 

Selection  of  Standard  Starting  Conditions .  As  was 
shown  in  the  previous  section,  it  takes  some  finite  time  for 
an  airfoil  at  angle  of  attack,  a,  suddenly  placed  into  motion 
to  build  to  a  steady-state  value  of  lift.  It  is  not  sur¬ 
prising,  then,  to  find  that  the  onset  of  constant  a  demon¬ 
strates  a  different  result  depending  on  the  time  delay  from 
onset  of  impulsive  motion  to  onset  of  constant  a.  The  dif¬ 
ferences,  however,  were  found  to  be  predictable,  and  thus 
separable,  as  the  following  will  show. 

To  determine  the  effect  the  initial  a  and  a  start  time 
have  on  the  vs.  a  curve  for  an  airfoil  at  constant  a,  a 

15%  thick  symmetric  Joukowski  airfoil  with  a  =  0.01  , 

★  .... 

At  =  0.02  was  started  at  various  initial  a's  and  allowed 

,  ★ 
to  build  lift  at  that  a  for  varying  lengths  of  time  t  .  For 

.  0  .  * 
initial  a  =  0  ,  one  can  see  from  Fig.  13  that  the  time  t 
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Effect  of  Start  Time  to  Begin  os 
15*  Joukowski  Airfoil,  a  =0?  a*=0.01 


at  which  the  a  begins  has  no  effect  on  the  vs.  a  curve. 

O 

For  initial  a  =  5  ,  the  vs.  a  curves  were  dependent 

★  • 

upon  the  value  of  t  at  which  the  a  was  begun.  However, 

as  can  be  seen  in  Fig.  14,  the  slope  of  the  vs.  a 

•  * 
curve  for  a  =  0.01  does  not  depend  upon  the  value  of  t 

•  ★ 

at  which  the  a  was  begun.  Choosing  At  =0.1  ,  the  15% 

thick  symmetric  Joukowski  airfoil  was  allowed  to  build  lift 

to  within  90%  of  steady-state  at  various  initial  a’s 

•  * 

before  starting  an  a  =  0.01  .  As  seen  in  Fig.  15,  the 

.  ,  o  o 

slope  of  the  vs.  a  curves  for  initial  a ’s  of  2  ,  4  and 
6°  are  all  approximately  equal.  The  dashed  lines  on  Fig.  15 
depict  the  vs.  a  curves  that  would  be  obtained  by  starting 
the  constant  -a  motion  at  full  steady-state  lift  values 
rather  than  the  90%  steady-state  lift  values  depicted  by 
the  solid  lines.  Note  that  the  initial  value  of  obtained 

o  o  o 

for  each  of  the  starting  angles  of  attack  of  2  ,  4  and  6 
is  the  same  amount  above  the  steady-state  curve,  and  is 
therefore  independent  of  initial  angle  of  attack.  This  ini¬ 
tial  value  of  will  be  called  the  ’jump'  condition.  Thus, 
by  the  foregoing  analysis,  vs.  a  curve  slope  effects  due 

to  the  vortex  wake  will  be  assumed  independent  of  initial  a 
★ 

and  t  . 

★ 

The  choice  of  At  also  shows  some  effect  on  the  vs. 

a  curve  and  was  investigated.  To  do  this,  a  15%  thick  sym- 

* 

metric  Joukowski  airfoil  at  a  =  0.01  was  run  at  At 
values  of  0.2,  0.1,  0.02  and  0.004.  Figure  16  depicts  a 

★ 

comparison  of  vs.  a  curves  for  these  four  values  of  At  . 
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0.02. 


* 

Note  that  the  slope  of  the  curves,  C,  ,  reduces  as  At 

a 

is  reduced,  but  the  reduction  is  negligible  below  At  = 

As  a  result  of  the  above  analysis  concerning  initial  a 

*  •  • 

and  t  at  which  a  is  begun,  all  constant-a  computer  runs 

assumed  an  initial  a  =  0°  and  t*  =  0  for  a  start-up.  A 
* 

standard  At  =0.02  was  chosen  as  a  reasonable  value  based 

upon  the  information  presented  in  Fig.  5  for  impulsive-start 

§  •  ★ 
motion  and  Fig. 16  for  constant-a  motion.  While  a  At  less 

than  0.02  would  produce  more  accurate  results,  the  increased 

★ 

computer  time  required  at  the  smaller  At  values  was  judged 
excessive  for  the  slight  increase  in  accuracy  that  could  be 
obtained . 

General  Effect  of  a  on  C,  .  To  determine  the  effect 
an  a  has  on  the  production  of  lift  on  an  airfoil,  a  15% 
thick  symmetric  Joukowski  airfoil  was  chosen  as  a  represen¬ 
tative  airfoil  shape.  Using  the  selected  values  of  initial 

O  *  • 

a  =  0  and  At  =0.02  ,  the  airfoil  was  subjected  to 

various  values  of  a*  ranging  from  0.005  to  0.035.  Figure  17 
depicts  the  vs.  a  curves  obtained  for  small  angles  of 

attack.  Comparing  with  the  vs.  a  curve  for  the  steady- 

•  ★ 

state  case,  one  can  see  that  as  the  value  of  a  is  increased, 
the  slope  of  the  Q„  vs.  a  curve,  C»  ,  is  reduced.  As  the 
motion  progresses  to  larger  values  of  a,  the  slopes  of  the 
curves  increase  slightly  (see  Fig.  18). 

Effect  of  Airfoil  Thickness  on  a  Effect .  The  general 
effect  a  has  on  the  production  of  lift  on  an  airfoil  has  been 
shown.  This  effect  was  shown  for  a  specific  airfoil  only. 


Fig.  18.  vs.  a  Slope  Change  as  a  Increases 
at  Constant  a 


To  determine  how  airfoil  thickness  may  influence  this  effect, 

several  symmetric  Joukowski  airfoils  of  varying  thickness 

were  subjected  to  the  same  a  conditions.  Once  again,  values 
*  o 

of  At  =  0.02  ,  initial  a  =  0  were  used.  Symmetric 

Joukowski  airfoils  of  7%,  15%  and  25.5 %  thickness,  as  well 
as  a  flat  plate  airfoil,  were  subjected  to  a  =  0.02.  A 
vs .  a  curve  can  be  plotted  for  each  of  these  airfoils. 
Plotting  the  average  slopes  of  these  curves,  C„  ,  versus  air- 
foil  thickness  ratio  t/c  (where  t  is  the  maximum  airfoil 
thickness),  one  can  determine  the  effect  of  airfoil  thickness 

on  the  vs.  a  curve  slope  reduction  due  to  a.  Figure  19 

.  *  . 

depicts  C .  vs.  t/c  for  a  =0.02  .  One  can  see  that  a 

has  a  greater  effect  on  lift  curve  slope  reduction  for  thin 

airfoils  than  for  thick  airfoils.  This  effect  is  consistent 

with  results  previously  presented.  Note  that  in  Fig.  9, 

where  a  =  0  ,  for  any  given  value  of  U^At/^c,  the  slope 

of  the  C ^/C %  curve  is  slightly  greater  for  the  25.5%  thick 

ss 

symmetric  Joukowski  airfoil  than  for  the  flat  plate.  Al¬ 
though  the  value  of  C^/C^,  is  less  for  the  airfoil  with 

ss 

thickness,  the  rate  at  which  C^/C,  is  increasing  is  greater. 

>s 

This  implies  that,  under  similar  a  conditions,  will  in¬ 
crease  at  a  faster  rate  for  a  thick  airfoil  than  for  a  thin 
airfoil.  Figure  19  confirms  that  conclusion. 

ft 

Effect  of  Airfoil  Camber  on  a  Effect .  In  much  the  same 
way  as  airfoil  thickness  effects  are  calculated,  airfoil 
camber  effects  can  also  be  explored.  Joukowski  airfoils  of 
15%  thickness  at  various  camber  ratios  were  subjected  to  an 
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a  =  0.02  .  As  before,  initial  a  was  0  ,  At  =  0.02 

Plotting  average  Cp  versus  camber  ratio  (maximum  camber/ 

*ot 

chord),  camber  effects  can  be  shown.  Figure  20  depicts  C 
vs.  camber  ratio  for  15%  thick  Joukowski  airfoils  of  various 
camber  ratios.  One  can  see  that  a  has  a  greater  effect  on 
lift  curve  slope  reduction  for  less  cambered  airfoils  than 
for  highly  cambered  airfoils. 


IV .  Conclus ion 


It  has  been  shown  that  as  an  airfoil  pitches  at  a  con¬ 
stant  a,  the  airfoil  trailing  vortex  wake  causes  the  slope 
of  the  vs.  a  curve  to  be  less  than  the  slope  of  the  vs.  a 
curve  for  steady-state  a.  The  greater  the  value  of  a  for  a 
given  airfoil,  the  greater  the  slope  reduction  of  the  C£  vs.  a 
curve  caused  by  the  vortex  wake.  This  effect  becomes  less 
pronounced  as  airfoil  thickness  increases.  Similarly,  the 
effect  is  also  less  pronounced  as  airfoil  camber  increases. 

Using  the  results  from  the  previous  section,  the  follow- 
ing  predictions  of  constant -a  effect  may  be  made. 

For  a  flat  plate,  the  reduction  in  C,  may  be  approx - 

a 

imately  calculated  by 

C  (a*)  =  — - 1 - +  2.2  (49) 

(a*  +  0 .00008) u 

(See  Fig.  21  for  a  comparison  of  this  prediction  with  numer¬ 
ical  data.)  This  prediction  may  be  approximately  corrected 
for  thickness  by  adding  a  correction  term  derived  from  Fig. 

19.  Thus 

C£  (a* , t/c)  =  (£)0,75  +  cz  (a*)  (50) 

a  a 

where  t/c  is  the  airfoil  thickness  to  chord  ratio  and 

C .  (a  )  is  the  C.  vs.  a  curve  slope  for  a  flat  plate  pre- 
a 

dieted  by  Eq.  (49).  A  further  approximate  correction  may 
be  made  for  camber  by  adding  another  correction  term,  de¬ 
rived  from  Fig.  20.  Thus 
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iction  of 


C„  (a  , t/c,m c/c)  =  e 


30(mc/c-0.09) 


+  C, 


(a  ,  t/c) 


(51) 


a 


where  mc/c  is  the  airfoil  camber  ratio  and  C„  (a  ,t/c)  is 
the  C ^  vs.  a  curve  slope  for  an  airfoil  with  thickness  pre¬ 
dicted  by  Eq.  (50). 

The  amount  that  increases  immediately  after  an  air¬ 
foil  begins  constant-a  motion  is  referred  to  as  the  •jump* 
condition.  This  value  for  a  flat  plate  can  be  predicted  by 


AC^ (a* )  =  3.47a*  (52) 

•  it 

where  AC^,(a  )  is  the  'jump'  condition  change.  Thickness 

*  * 

effects  on  AC^(a  )  can  be  approximated  by  the  equation 

AC^ (a* , t/c)  =  [1  +  2i£Z£l]ACi(d*)  (53) 

*  *  . 

where  t/c  is  the  airfoil  thickness  ratio  and  AC^(a  )  is  the 
'jump'  condition  for  a  flat  plate  defined  by  Eq .  (52).  A 
final  approximate  correction  to  the  ’jump'  condition  can  be 
made  by 

•  ★  • ★  mo 

AC^(a  t/c,  mc/c)  =  AC^(a  ,  t/c)  -  1.3(  — )  (54) 

where  mc/c  is  the  airfoil  camber  ratio  and  AC^,(a  ,  t/c)  is 
as  defined  by  Eq .  (53). 
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V .  Recommendation 


The  assumption  that  the  trailing  vortex  wake  of  an  air¬ 
foil  undergoing  a  constant  rate  of  change  of  angle  of  attack 
has  a  negligible  effect  on  the  production  of  lift  on  the  air¬ 
foil  is  not,  in  general,  valid.  Although  the  effect  is  not 
large  (see  Eq.  (49)),  it  should  be  accounted  for  in  the  in¬ 
vestigation  of  dynamic  stall  of  airfoils.  The  methods  de¬ 
veloped  by  Docken  (4)  and  Lawrence  (5)  could  be  modified  to 
include  the  techniques  presented  in  this  thesis  to  more  ac¬ 
curately  predict  the  pbtential  flow  field  about  a  pitching 
airfoil  at  any  instant  in  time.  Incorporating  the  calcula¬ 
tion  of  wake  vortex  effects  outlined  in  this  thesis  into 
Lawrence's  work  would  significantly  contribute  to  the  solu¬ 
tion  of  the  dynamic  stall  problem  for  an  airfoil  undergoing 
a  constant  rate  of  change  of  angle  of  attack. 


I 
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APPENDIX;  Computer  Program 


This  program  computes  circulation,  pressure  difference 
distribution,  vorticity  d istribution ,  coefficient  of 
lift,  and  trailing  vortex  wake  shope  for  a  2-D  JouKowsKi 
airfoil  in  an  incompressible,  inviscid  free  stream  at 
angle  of  attack.  The  angle  of  attack  may  be  a  constant 
value,  or  it  may  be  changed  at  a  constant  rate  for  the 
number  of  time  steps  desired.  All  output  values  are 
computed  assuming  a  trailing  vortex  wake  made  up  of  dis¬ 
crete  point  vortices  of  constant  strength,  each  of  which 
influences  the  motion  of  all  the  other  vortices  and  the 
flow  about  the  airfoil.  For  the  constant  rate-of-chanae 
of  angle-of-attack  case,  coefficient  of  lift  can  be 


found  as  a  function  of  the  rate  of  change  of  angle  of 


attack  . 

Variables  in  the  program  are  defined  as  follows 

a  1  f  a 

- 

angle  of  attack 

alfaO 

- 

initial  angle  of  attack 

alf dot 

- 

time  rate  of  change  of  angle  of  attack 

beta 

- 

the  beta  parameter  of  a  Joukowski  airfoil 

calf  a 

- 

COSINE  of  alfa  • 

chord 

- 

airfoil  chord  length 

cl 

- 

coefficient  of  lift 

c  iss 

- 

steady-state  coefficient  of  lift 

count! 

an  integer  counter  used  to  determine  which 
time  steps  will  record  output  in  certain 
f  i .  e  s 

c theta 

- 

COSINE  of  theta 

d  a  1  f  a 

- 

incremental  change  of  alfa 

CiOiCD 

array  of  incremental  values  of  coefficient 
of  pressure  along  the  airfoil  chord 

deid 

distance  on  the  x-axis  in  the  cylinder 
plane  behind  the  cylinder  where  the  first- 
shed  vortex  is  placed 

del gam 

array  of  values  of  strengths  of  gamma  for 
each  individual  vortex  pair 

daams 

— 

a  sum  of  vortex  strengths 

dsl 

incremental  distance  along  airfoil  lower 
surface 

dsu 

" 

incremental  distance  along  airfoil  upper 
surface 

dt 

- 

incremental  unit  of  time 

DZD 

' 

complex  number?  derivative  of  the 

Joukowski  transformation 

dccrc 

~ 

magnitude  of  DZD 

e  t.  a 

x- value  of  trailing  vortex  position 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

r 

w 

c 

c 

c 

c 

c 

c 

c 

w 

c 

c 

c 

c 

c 


c 

n 

w 

c 

/*■» 

u 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


eta2 
gamma 
i  t  .j  *  k  *  1 
intg  rl 


lastdt 
mag  2 

malfa 

maxdt 

maxt 

MU 


pi 

RHO 

RHOP 

sal  fa 
ssaam 
stheta 
sumsq r 

sumu  * sumv 


t 

theta 
tOtU  r totV 


IJ  *  V 

ua  *  va 


usurf * vsurf 

utheta 
vord is 
vsqare 


Xr  y 

xvort . yvort 

seta 

r«tc2 


-  x-vaiue  of  trailing  vortex  image  position 

-  airfoil  circulation 

-  integer  values  used  in  iterations 

-  value  of  the  integration  of  velocity 
differences  between  upper  and  lower 
airfoil  surfaces 

-  last  time  increment  at  which  alfa  changes 

-  distance  from  a  vortex  to  the  center  of  the 
cylinder  in  the  cylinder-centered  plane 

-  maximum  volue  of  alfa 

-  first  time  increment  at  which  alfa  changes 
from  alfaO 

-  last  time  step 

-  complex  number?  distance  between  the  origin 
in  the  d isp laced-cy 1 inder  plane  and  the 
center  of  the  cylinder 

-  the  constant  3.14159265 

-  complex  number?  a  position  in  the  cylindei — 
centered  plane 

-  complex  number;  a  position  in  the  displaced- 
cylinder  plane 

-SINE  of  alfa 

-  steady  state  value  of  circulation 

-  SINE  of  theta 

-  the  square  of  the  distance  of  a  trailing 
vortex  from  the  origin  in  the  cylinder  plane 

-  sum  of  the  velocities  on  the  cylinder  in 
the  x  and  y  directions*  respectively *  in¬ 
duced  by  trailing  vortices  and  their  images 

-  integer  counter  for  number  of  time  steps 

-  angle  measured  counterc locKwise  from  the 
x-axis  in  the  cylinder  plane 

-  sum  of  velocities  at  a  vortex  location  in 
the  x  and  y  directions*  respectively t  in¬ 
duced  by  trailing  vortices  and  their  images 

-  velocities  at  a  vortex  location  in  the 
x  and  y  directions*  respectively 

-  velocities  at  a  vortex  location  in  the 

x  and  y  directions*  respectively*  in  the 
airfoil  plane 

-  velocities  on  the  cylinder  in  the  x  and  y 
directions*  respectively 

-  velocity  tangent  to  the  cylinder 

-  vorticity 

-  velocity  on  upper  surface  of  airfoil 
squared  minus  velocity  on  lower  surface 
of  airfoil  squared 

-  position  on  airfoil 

-  position  of  a  trailing  vortex 

-  complex  number?  a  position  on  the  airfoil 

-  y  value  of  trailing  vortex  position 

-  y  value  of  trailing  vortex  image  position 
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zt  -  distance  along  x-axis  frori  origin  to 

cylinder  in  displaced -cylinder  plane 
ZZ  -  complex  number;  a  position  in  the  wake  in 

the  airfoil  plane 

files; 


INPUT 

OUTPUT 


PRESD 

VORT 

WAKE 


unformatted  list  of  input  variables 

list  of  Cl  vs.  Cl  steady-state  for  angle  of 

attack  and  time  step 

pressure  distribution  at  specified  time 
vorticity  distribution  at  specified  time 
position  of  trailing  vortices  at  specified  time 


DIMENSION  seta  <  201 , 201 ) , eta < 201 ,201 ) , etu2 ( 20] ,201) , 
fzetaz < 201,201 > ,u (201, 201) ,v(201 ,201 ) ,delgam(201 )  , 

+x(-l 80 J1 80) ,y (-180J 180) , utheta < -180 J 180) ,uo (201,201 ) , 
+va ( 201 ,201 ) ,xvort(201 ,201 > ,yvort (201 ,201 ) , 
rvordis<0: 176) ,delcp<0J 176) , intgrl (0 J201 ,0 J 180) 

INTEGER  i  ,  .j  ,  maxt ,  t  ,k  ,  1 ,  maxd  t  ,countt 

REAL  delgam ,deld , pi »c theta , stheta ,ssgam,dsu,dsl,delcp, 
fusurf , vsurf , theta, u theta, gamma, malf a, al fa, sumu,sumv, 
fdalf a,dgams,calfa,salfa,sumsqr,alfdot,dt,dzdro,ua,va, 
fbeta , zt ,alfa0,clss,cl , x , y , vsqare, xvort ,yvort , lastdt , 
+et a2 , zeta , zeta2 , u, v, chord , intg  rl , totu , to tv, vordis, 
+mag2» eta 

COMPLEX  MU , RH0P , Z , ZZ , RH0 , DZP 
OPEN  < 15, FILE*' INPUT' ) 

REWIND  15 

OPEN  (16,  FILE* 'OUTPUT') 

REWIND  16 

OPEN  (17. FILE* 'PRESD') 

REWIND  17 

OPEN  ( 18,FILE='WAKE' ) 

REWIND  13 

OPEN  ( 1?,FILE='V0RT' ) 

REWIND  1? 
pi=3. 14159265 
C 

C  Initialize  delta  alpha,  delta  d,  max  t.  Compute  steady 
C  state  gamma. 

C 

10  CONTINUE 

READ( 15,#,END=400)  beta,alfaO,dalfa,deld , maxd t, maxt, 
+lastdt,zt,alfdot,dt 
malf  a=dalf a#< lastdt-maxdt ) +alf aO 
ssgam=4*pi#SIN ( malf a#p i/ 180 ) 

WRI TE( 16, 70 )  zt , beta , da  If a , ceid , ssgan . alf  dc t , dt 
WRITE ( 17, 72 )  zt ,beta , dal f a , aeid , ssgam, alf dot ,dt 
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WRITE ( 18, 72)  zt ,beta , dalfa,deld , s£gam,alf dot , dt 

WRITE ( 1 9 ,72 )  zt ,  beta ,dalfa,deld , ssgam , alf dot » dt 

beta=beta*pi/180 

da  If a=dalf a*p i/180 

alfa0=alfa0*pi/180 

MU=CMF'LX<zt-COS(beta) ,SIN<beta>  > 

Calculate  coordinates  of  points  on  the  airfoil,  (x,y) 

DO  15  i=-180, 180 
theta=i#pi/180 

RH0P=CMPLX<C0S<theta-beta> ,SIN< theta-beta) >+MU 
Z=RH0P+zt**2/RH0P 
x < i ) =REAL <  Z ) 
y ( i )=AIMAG(Z) 

15  CONTINUE 

chord =x  ( 0 )  -x  ( -1 80+2*beta ) 

DO  12  i=l,maxt 
del gam  (  i  )  =0 . 0 

12  CONTINUE 

BO  13  i=176,0,-4 
intgrl (0, i )=0, 

13  CONTINUE 

Begin  stepping  in  time,  inserting  a  new  vortex  pair  at 
each  time  step. 

coun tt=0 
BO  300  t=l,maxt 
countt=countt+l 
intgrl <t,l 80 )=0. 
daams=0 .0 
cl=0. 

IF  ( t .GE.maxdt .and  » t .LE. iastdt )  THEN 
a I f o= < t-maxdt ) #dal fatal f aO 
ELSE  IF(t.LT.inoxdt)  THEN 
alf a=alf aO 
END  IF 

calf a=C0S ( alf a ) 
salf a=SIN<alfa) 

Insert  new  vortex  pair,  and  update  position  of  all  other 
vortex  images, 

zeta (t,l)=*deld+l. 
eta  < t , 1 )=0» 

xvort ( t , 1 ) =chord#dt*COS ( 2#beta> /2+x (0) 
yvort<t,l)=chord*dt*SIN(2*beta>/(-2> 

DO  20  , j  =  1 , t 

zeta2 ( t  ,.j )  =zeta  ( t ,  J  )  /  <  zeta  ( t ,  J  )  *#2+eta  ( t ,  ,j )  **2 ) 
eta2  (  \  ,  .J )  =eta  <  t  *  .j  )  /  ( zeta  ( t ,  .j )  *#2reta  <  t ,  J  )  ##2  > 

20  CONTINUE 
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Calculate  strength  of  newly  shed  vortex  at  time  t. 

This  is  done  by  satisfying  the  Kutta  Condition  at  the 
trailing  edge  while  Keeping  total  circulation  in  the 
field  equal  to  zero. 

delgam  < 1 ) =4*pi*salf a*< zeta < t , 1 > **2+eta ( t , 1 > **2+1-2* 

+  zeta(t, 1) )/<zeta(t,l)**2+eta<t,l>**2-l> 

IF(t.GE.2>  THEN 
DO  21  i=t,2,-l 

dgams=dganis+delgain(i  >*<zeta<t,  i  >**2+eta(  t,  i  >**2-1 ) 
+  / ( zeta<  t ? i>**2+eta<tfi) **2+l-2*zeta<  t ,  i )  >* 

+  delgamd  )/(4*pi*salfa) 

21  CONTINUE 
END  IF 

delgam< 1 )=delgam( 1 )-dgams 

Compute  velocities  and  circulation  around  the  cylinder 
at  time  t. 

DO  200  i=-179, 178,2 
theta=i*pi/180 
ctheta=COS ( theta ) 
stheta-SIN< theta) 
sumu=0 . 0 
sumv=0 . 0 
DO  30  k=t, 1 ,-l 

sumu=delgam(k )/2/pi*( <eta(t*k )-stheta)/C (ctheta- 
+  zeta ( t . K ) )  **2+  <  s theta-eta ( t  f  k  ) ) **2  >  +  <  stheta- 

+  eta2 ( t f K ) )/( (ctheta-zeta2(t,K ) >**2+(stheta- 

+  eta2 < t fk ) )**2 ) ) +sumu 

sumv-de] gam(K ) /2/pi*< (ctheta-zeta ( t » K  > )/ < ( c theta - 
+  zeta ( t » K ) ) **2+ <  stheta-eta ( t , k ) )**2)-( (ctheta- 

+  zeta2(t , K ) )/(  <ctheta-zsta2< t ,K > ) **2+ ( stheta- 

■f  sta2(tfk ) >**2> ) )+sumv 

30  CONTINUE 

usurf =2* (calf  a*stheta**2-c theta*stheta*salfa )+sumu 
vsurf =2* ( sal f a*ctheta**2-ctheta*stheta*calfa ) +sumv 
utheta( i )=SGRT  <usurf**2+vsurf**2) 

IF(i.NE.O)  THEN 

RHQP=CMPLX ( COS ( theta-beta ) *  SIN ( theta-beta ) >+MU 
DZD=l-zt**2/RH0P*#2 
dzdro*6BS(DZD) 
utheta< i)=utheta( i )/dzdro 
END  IF 

IF( < theta- <alfa+2*beta) ).LT«-l*pi)  THEN 
u theta  < i )=-l*utheta ( i > 

END  IF 

200  CONTINUE 

gamma=0 . 

DO  210  i  =  1 r  t 

gamma^gaiTiiTiatdelcamC  i  ) 
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210  CONTINUE 


Compute  velocity  at  each  discrete  vortex  location  due 
to  all  other  vortices  and  free  stream. 

DO  50  1=1, t 
totu=0.0 
totv=0»0 
DO  40  K=l,t 

IF(k.EQ.l)  THEN 

totu=delgam<K )/2/pi#< <  eta( t ,K  >-eta2< t,k ) )/( ( 

+  zeta ( t , k  > -zeta2 <  t ,  k  > ) **2+ ( eta  <  t , k ) -eta2 ( t , k  > ) 

+  **2))+totu 

totv=delgam<k )/2/pi#( <zeta2(t,k  >-zeta< t,k ) )/( ( 

+  zeta ( t , k ) -zeta2 <  t ,  k  >  > **2+  <  eta ( t , k ) -eta2 ( t , k ) ) 

+  *#2)>-i-totv 

ELSE 

totu=delgam<K )/2/pi# ( ( eta< t , 1 >-eta2( t ,K ) )/< ( 

+  zeta ( t , 1 ) -zeta2 ( t , k ) ) **2+ ( eta ( t , 1 ) -eta2 ( t , K  > ) 

+  **2 )  +  ( eta  <  t , k ) -eta  <  t , 1 ) > / ( < zeta ( t , 1 > -zeta ( t , K  > ) 

+  **2+<eta<t,l)-eta(t,k>  >**2) >+totu 

totv=delgam<k)/2/pi*<  <zeta2( t ,k )-zeta( t, 1 ) )/( ( 

+  zeta ( t , 1 ) -zeta2 <  t , k  > ) **2+ ( eta ( t ,  1  > -eta2 ( t , k ) ) 

+  ##2 )+ (zeta< t , 1 )-zeta  <  t ,k ) >/ ( (zeta ( t , 1 >-zeta ( t ,k ) 

+  )**2+<eta(t,l)-eta(t,k>  >**2) >+totv 

END  IF 

40  CONTINUE 

sumsq  r'-zeta  ( t ,  1 )  **2+eta  ( t ,  1  >  *#2 

u<t,l)=calfa*<l+<eta(t,l>**2-zeta(t,l)**2)/surosqr 
+  *#2)-2#zeta  ( t ,  1  >#eta  <  t ,  1 )  *salfa/sumsqr**'2+totu 

v<  t, 1 )=salfa#< l  +  (zeta( t, 1 )  #*2-eta ( t , 1 >##2>/sumsqr 
r  *#2 )-2#zeta  <  t ,  1  >#eta  <  t ,  1 )  *calf  a/sumsqr*)|c2+totv 

50  CONTINUE 

Calculate  pressure  distribution  and  unsteady  aerodynamic 
force  on  the  airfoil. 

DO  32  i=176 , 0 , -4 
theta  =i#ai/180 
ctheta=C0S< theta ) 
stheta=SIN( theta) 

d  su*S0RT ( ( x ( i ) -x ( i+4 ) ) **2+  <  y ( i > -y  < i +4 ) ) **2 ) *2/chord 
dsl“SGRT  < (x(-l*i)-x<-l*(i+4)))**2+<y(-l*i>-y(-l# 

+  <i+4> ) >**2)*2/chord 

intgrl  ( t ,  i )  =utheta  (  i+2)#dsu-utheta<  -1*<  i+2)  )*dsl 
+  tintgrl <t, i+4) 

vsqare=utheta<  i  )*#2-utheta(-l#i  )  *#2 
IF(t.EG.l)  THEN 

delcp  < i )=vsqare+2*intgrl (1 ,i >/dt 
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ELSE  IF(t.EQ.2>  THEN 

delcp  (  i  )=vsqare+2*<  intgrl  <2,  i  )-intgrl  <  1 ,  i  )  )/dt 
ELSE  IF(t.EQ.moxdt)  THEN 

delcp ( i )=vsqare+2* intgrl (maxdt, i >/dt 
ELSE  IF ( t ♦  EQ . maxdt+1 )  THEN 

delcp ( i )=vsqore+2*( intgrl (maxdt+1 , i )- 
intgrl (maxdt, i ) )/dt 

ELSE 

delcp ( i >=vsqare+< intgrl ( t-2, i )+3*intgrl ( t , i )-4* 

+  intg rl < t-1 , i ) )/dt 

END  IF 

IF ( i . EQ ♦ 176 >  THEN 

cl=cl+delcp(176)*(x(174)-x(180) ) /chord 
ELSE  IF(i.EG.O)  THEN 

c l=c 1+delcp ( 0 ) # (x ( 0) -x ( 2) ) /chord 
ELSE 

cl=c 1+delcp ( i )*(x ( i-2) -x( i+2) ) /chord 
END  IF 

22  CONTINUE 

Calculate  Cl. 

c lss=8#p i/chord#salf a 
IF  <  t . EQ . 3 )  THEN 

WRITE ( 17  * ' ( 'Delta  Cp ,  ,fI3f*  Vortices,  t=,.F7.4)'> 

+  t,t*dt 

WRITE ( 19 ,'(* Aortic ity  Distribution  *,I3» 

+  *  Vortices,  t=’,F7.4)')  t,t*dt 

DO  34  i=176 , 0 , -4 
vordis( i )=utheta( i )-u theta ( -l#i ) 

WRITE ( 19, 90 )  <  x< i )+chord-2*zt ) *2/chord , vord is ( i ) 
WRITE ( 17,90)  ( x ( i ) +chord-2*zt )#2/chord ,-l#delcp(i) 

34  CONTINUE 
END  IF 

IF ( c oun 1 1/ 1 0 . GT . 0  >  THEM 

WRITE(17, '( 'Delta  Cp,*,I3,*  Vortices,  t=*,F7.4)') 

+  t,t*at 

WRITE(19, ' ( 'Vorticity  Distribution  ’,13, 

+  *  Vortices,  t=*»F7.4>')  t,t#dt 

DO  35  i=176 , 0 , -4 
vord is < i)=utheta( i ) -utheta( -l#i ) 

WRITE ( 19, 90)  (x< i ) +chord-2#zt )#2/chord , vordis( i ) 
WRITE( 17,90)  <x< i >+chord-2#zt )#2/chord ,-l*delcp  < i ) 

35  CONTINUE 
countt=0 

END  IF 

WRITE < 16,80)  t, (a  If a-b eta )*1 80/pi , (xvort(t , t)-x(O) ) 

+  #2/chord ,yvort(t,t)#2/chord , gamma, clss, cl 

Hove  each  vortex  to  its  new  location  in  the  flow  field. 
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DC  60  R=t ,1,-1 

mag 2=SQRT( seta ( t  r  k  )*#2+eta  ( t ,  k  >##2) 

RHQPs=mag2*CMPLX  ( COS  ( ACOS  <  seta  ( t  ,  K  ) /mag2  >  -beta )  , 

SIN  ( ASIN  (  eta ( t  ,  K ) /mag2  > -beta ) ) -f  MU 
DZD=l-st**2/RHQP**2 
d2dro=ABS(DZD> 
ua ( t ,  K )=u ( t ,k )/dzd ro 
va(t,k )=v(t,k >/dzdro 

x vo  rt  (  t +1 , K +1  >  =xvo  rt  ( t ,  K )  +ua  <  t  ,  K  >  #cho  rd*d  t/2 
yvort (t+1 ,K  +  l)=yvort (t,K  >+va(t»K ) #chcrd*dt/2 
ZZ=CMPLX(xvort(t+l,K+l> ,yvort(t  +  l,K+l>  > 

RHO=< (ZZ+  SORT ( ZZ**2-4*zt**2 )  ) /2-MU)* 

+  CMPLX(pOS(beta) ,SIN(beta> ) 

seta (t+1 ,  K  +  l )=REAL(RHO) 
eta(t+l,K+l)=AIMAG(RHO) 
del gam (K+l ) -del gam (K ) 

60  CONTINUE 

IF (t .EQ.5)  THEN 
GO  TO  61 

ELSE  IF  (t.EQ.maxdt)  THEN 
GO  TO  61 

ELSE  IF < t *EB. lastdt )  THEN 
GO  TO  61 

ELSE  IF (t.EG.maxt)  THEN 
GO  TO  61 

ELSE 

GO  TO  300 

END  IF 

61  WRITE ( 18, ' ( 'Wake  Vortex  Locations  from  Trailing*, 

+*  Edge  (1/2  c  =  1)*/*  Xvort * , 5X , * Yvort ’ ) ' > 

DO  65  i  =  l,t 

WRITE ( 18, 95 )  ( xvort <  t » i )-x( 0) ) #2/chord , 

+  yvort ( t , i > *2/chord 

65  CONTINUE 

00  CONTINUE 

70  FORMAT (///,' AIRFO IL  DATA  :',//,'Zeta  trailing  edge  :', 
+F7 .  4 » '  Beta  (deg rees ) 5 ' ,F6 . 3 DYNAMIC  PARAMETERS ' , 
+//, 'Delta  Alpha  (degrees) {' ,F6.3, '  Delta  Vortex  ", 

+  *  Distance  *'»  F6 . 4  r  /,'  Steady  State  Gamma  t  '»F7.5, 

+  '  Alpha  Dot J',F8»5,'  Delta  Time: ' ,F5.3,//, 'Time' , 10X, 
t'Starting  Vortex' , 17X, 'Cl ',/, 'Step  Alpha  X  ', 

+  '  Y  Gamma  Steady  State  CIS/) 

72  FORMATS'//,  'AIRFOIL  DATA  :',//,'Zeta  trailing  edge 

+F7 ♦  4 ,  '  Beta  (degrees) J ',F6. 3,//, 'DYNAMIC  PARAMETERS', 
+//, 'Delta  Alpha  (degrees) J ',F6. 3, '  Delta  Vortex  ', 

t'Distance: ' ,F6»4,/, 'Steady  State  Gamma:  SF7.5, 

+'  Alpha  Dot* ' ,F8»5, '  Delta  Time: ' ,F5.3,/) 

80  FORMAT ( 13 , 3X ,F6 *3, 2X,F7 »4,F8.4,1X,F8»5,5X,FB*5,2X,FB»5) 
90  FORMAT (FS»3,4X,F10*5) 

95  FCRMAT(F8.4,F11.4) 

00  END 
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